% Program to estimate Swank-Steinmo Reanalysis for JTP article
% Created by Jude Hays November 22, 2006

clear all;
A=wk1read('C:\Documents and Settings\Administrator\Desktop\JTP Volume\swank_reanalysis\TAXDATA4.wk1',1,0);
W1=wk1read('C:\Documents and Settings\Administrator\Desktop\JTP Volume\swank_reanalysis\TDW5.wk1');
%M=open('C:\Documents and Settings\jchays\Desktop\Research\Book\chapters\Chapter 2\Chapter_2_Data\M_GC.mat');
%W_LR_1=wk1read('C:\Documents and Settings\jchays\Desktop\Research\Book\chapters\Chapter 2\Chapter_2_Data\W_GC_LR.wk1');
nobs=195;
% normalize W if still necessary
W_LR=normw(W1);
I_T = eye(15);
W1_TSCS = kron(I_T,W1);
W=(normw(W1_TSCS)); %normw
y=A(:,[3]); % column number in the data matrix that corresponds to the dependent variable
x=A(:,[4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33]);%,36,37,38,39,40]); % column numbers in the data matrix that correspond to the independent variables 8,9,10,11,12,
xconstant=ones(nobs,1);
info.rmin=-.9;
info.rmax=.9;
info.lflag=0; % 0 required for exact results
results=sar(y,[xconstant x],W,info);
vnames=strvcat('y','constant','l1_dv','l1_capcon','l1_trade','prof','dominv','l1_ltunemp','psdebt','l1_old','l1_grow','l1_leftc','l1_cdemc','c2','c3','c4','c5','c6','c7','c8','c9','c10','c11','c12','c13','y823','y845','y867','y889','y901','y923','y945');,%'unemp','infl'
prt(results,vnames);
B_TEMP=[results.beta;results.rho];
BSTD_TEMP=[results.bstd;results.pstd];
VCVM_TEMP=[results.vcvm];
VCVM = [VCVM_TEMP(32,32) VCVM_TEMP(32,2) VCVM_TEMP(32,5);  
VCVM_TEMP(2,32) VCVM_TEMP(2,2) VCVM_TEMP(2,5);
VCVM_TEMP(32,5) VCVM_TEMP(2,5) VCVM_TEMP(5,5)];

VCVM2 = [VCVM_TEMP(32,32) VCVM_TEMP(32,5);  
VCVM_TEMP(32,5) VCVM_TEMP(5,5)];

%Long-Run Steady-State Spatio-Temporal Effects
i_n=eye(13);
estcf1_LR = i_n - B_TEMP(32,1)*W_LR - B_TEMP(2,1)*i_n;
estcf_LR = (B_TEMP(5,1))*inv(estcf1_LR); %conditionial effect of LT_Unemployment
d_cf_rho1_LR = estcf_LR*W_LR*inv(estcf1_LR);
d_cf_rho_LR = d_cf_rho1_LR(1:13,13); 
d_cf_phi1_LR = estcf_LR*inv(estcf1_LR);
d_cf_phi_LR = d_cf_phi1_LR(1:13,13);
d_cf_beta1_LR = inv(estcf1_LR);
d_cf_beta_LR = d_cf_beta1_LR(1:13,13);
d_cf_TEMP_LR = [d_cf_rho_LR d_cf_phi_LR d_cf_beta_LR];
d_cf_LR = d_cf_TEMP_LR;
delta1_LR = d_cf_LR*VCVM*d_cf_LR';
delta_LR = diag(delta1_LR);
cf_LR = estcf_LR;

%Short-Run Spatial Effects
i_n=eye(13);
estcf1_SR = i_n - B_TEMP(32,1)*W_LR;
estcf_SR = (B_TEMP(5,1))*inv(estcf1_SR); %conditionial effect of LT_Unemployment
d_cf_rho1_SR = estcf_SR*W_LR*inv(estcf1_SR);
d_cf_rho_SR = d_cf_rho1_SR(1:13,13); 
d_cf_beta1_SR = inv(estcf1_SR);
d_cf_beta_SR = d_cf_beta1_SR(1:13,13);
d_cf_TEMP_SR = [d_cf_rho_SR d_cf_beta_SR];
d_cf_SR = d_cf_TEMP_SR;
delta1_SR = d_cf_SR*VCVM2*d_cf_SR';
delta_SR = diag(delta1_SR);
cf_SR = estcf_SR;

%Response Functions
nobs = 13*15;
    m1 = eye(nobs-13,nobs-13);
    m2 = zeros(13,nobs-13);
    m3 = zeros(nobs,13);
    m4 = [m2 ; m1];
    M = [m4 m3];
i_nt=eye(nobs);
estcf1_RF = i_nt - B_TEMP(32,1)*W - B_TEMP(2,1)*M;
estcf_RF = (B_TEMP(5,1))*inv(estcf1_RF); %conditionial effect of LT_Unemployment
d_cf_rho1_RF = estcf_RF*W*inv(estcf1_RF);
d_cf_rho_RF = d_cf_rho1_RF(1:nobs,6); 
d_cf_phi1_RF = estcf_RF*inv(estcf1_RF);
d_cf_phi_RF = d_cf_phi1_RF(1:nobs,6);
d_cf_beta1_RF = inv(estcf1_RF);
d_cf_beta_RF = d_cf_beta1_RF(1:nobs,6);
d_cf_TEMP_RF = [d_cf_rho_RF d_cf_phi_RF d_cf_beta_RF];
d_cf_RF = d_cf_TEMP_RF;
delta1_RF = d_cf_RF*VCVM*d_cf_RF';
delta_RF = diag(delta1_RF);
cf_RF = estcf_RF;